#include "SphericalHarmonics.h"

template<typename T>
const T SHFunction<T>::PI = static_cast<T>(std::acos(-1.0));

template<typename T>
SHFunction<T>::SHFunction(int bands)
{
	nbands = bands;
	values.assign(bands * bands, T(0));
}

template<typename T>
T SHFunction<T>::random()
{
	return static_cast<T>(rand()) / static_cast<T>(RAND_MAX);
}

template<typename T>
void SHFunction<T>::random(T &theta, T &phi)
{
	T x1 = random();
	T x2 = random();
	theta = 2 * std::acos(sqrt(x1));
	phi = 2 * PI * x2;
}

template<typename T>
T SHFunction<T>::realY(int l, int m, T theta, T phi)
{
	if(m >= 0)
		return boost::math::spherical_harmonic_r(l, m, theta, phi);
	else
		return boost::math::spherical_harmonic_i(l, m, theta, phi);
}

template<typename T>
T SHFunction<T>::evaluate(T theta, T phi) const
{
	T value = T(0);
	for(int l = 0; l < nbands; ++l) for(int m = -l; m <= l; ++m)
		value += values[l * l + l + m] * realY(l, m, theta, phi);
	return value;
}

template class SHFunction<float>;
template class SHFunction<double>;
